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Although the hypothesis that nestedness determines mutualistic ecosystem dynamics is accepted in general, 
results of some recent data analyses and theoretical studies have begun to cast doubt on the impact of 
nestedness on ecosystem stability. However, definite conclusions have not yet been reached because previous 
studies are mainly based on numerical simulations. Therefore, we reveal a mathematical architecture in the 
relationship between ecological mutualistic networks and local stability based on spectral graph analysis. In 
particular, we propose a theoretical method for estimating the dominant eigenvalue (i.e., spectral radius) of 
quantitative (or weighted) bipartite networks by extending spectral graph theory, and provide a theoretical 
prediction that the heterogeneity of node degrees and link weights primarily determines the local stability; 
on the other hand, nestedness additionally affects it. Numerical simulations demonstrate the validity of our 
theory and prediction. This study emphasizes the importance of ecological network heterogeneity in 
ecosystem dynamics, and it enhances our understanding of structure-stability relationships. 



Ecological communities consist of a number of species that are connected via interspecific interactions, such 
as trophic and mutualistic relationships. Their structure and dynamics are significant in ecology, because 
they are not only important in the context of basic scientific research but also in the context of applied 
ecology'"^. In particular, the nature of the structure-stability relationship is a long-standing question'"^. In this 
context, eigenvalue analysis, which reflects how a dynamical system resting at equilibrium responds to perturba- 
tions, has been widely used when discussing the local stability of ecosystems' '. 

Ecological communities are often represented as networks' **. Previous network analytical studies have revealed 
that plant-animal mutualistic systems exhibit a non-random structural pattern: nested architecture', which 
occurs when the interaction pairs of a certain (specialist) species form a subset of those of another (generalist) 
species in a hierarchical fashion. Importantly, nestedness is believed to influence ecological dynamics. Nestedness 
may minimize competition and increase biodiversity in mutualistic networks', and it emerges as a result of an 
optimization principle aimed at maximizing species abundance'". Moreover, nestedness promotes the resilience 
and persistence of mutualistic networks, although it inhibits them in predator-prey networks". The extinction of 
species, which contributes to the emergence of nestedness in mutualistic networks, decreases the persistence of 
ecosystems'^. 

However, several recent studies have begun to cast doubt on the importance of nestedness. For example, nested 
architecture can be more easily acquired than previously thought'^ Staniczenko et al.' demonstrated that nest- 
edness can be concluded in binary networks (i.e., presence 1, or absence 0, of a given link), but not in quantitative 
or weighted networks. Moreover, they used a spectral graph approach to show that the largest eigenvalue of a 
community matrix determines nestedness. James et al.'* reported that biodiversity in mutualistic communities is 
described by the number of mutualistic partners a species has (i.e., node degree) rather than nestedness. Jonhson 
et al.'^ mathematically showed that nestedness can result from the heterogeneity of degree distributions. 

Heterogeneities of degree distributions for binary networks and strength (i.e., the sum of neighbors' link 
weights) distributions for quantitative networks are significant properties in wide-ranging real-world networks, 
including mutualistic networks"", and they influence the dynamics on the networks such as synchronization"'"*, 
controllability", epidemics of infectious diseases"^", and so on (see^' for a review). 

However, it is poorly understood whether nestedness or the heterogeneity of degree (or strength) distributions 
is most important for explaining the dynamics of mutualistic ecosystems, because the previous studies reviewed 
above were mainly based on numerical simulations. Therefore, in this study, we investigate the impact of 
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structural properties (e.g., heterogeneity and nestedness) on local 
stability (i.e., the dominant eigenvalue) of mutualistic ecological 
communities in greater detail. In particular, we propose a theoretical 
method — the Semicircle Plus Twins Method — for estimating the 
dominant or largest eigenvalue of both binary and quantitative 
bipartite networks. We do this by extending the superelliptical law 
for large-scale unipartite binary networks"'^, because it is not directly 
applicable to either binary or quantitative bipartite networks. We 
also confirm the validity of this theoretical method using numerical 
simulations. In addition, we apply this estimation method to 40 
empirical mutualistic networks, and evaluate the effect of structural 
properties on the local stability of mutualistic communities using 
nuU model analysis. The proposed method and numerical simula- 
tions show that the local stability of mutualistic ecosystems is deter- 
mined by heterogeneous properties rather than topological 
nestedness such as NODF and WNODF (see Methods for details). 

Results 

Theoretical estimation of the dominant eigenvalue of bipartite 
networks. We here introduce analytical methods for estimating 
the dominant eigenvalue of both binary and quantitative bipartite 
networks and reveal a relationship between the dominant eigenvalue 
and structural properties. 

Thus far, several important laws for estimating eigenvalues have 
been reported. In particular, Wigner's semicircle law^^'^'* is well 
known. 

Wigner's semicircle law states that the density of real eigenvalues, 
obtained from a large symmetrical random matrix A, follows a semi- 
circular (half-moon) distribution when A satisfies the following con- 
straints: 1) the elements are sampled from a zero mean distribution 
(i.e., E[Aij) =0) and 2) A is dense (i.e., as s ^ ^, sc ^ o=, where s 
corresponds to the number of nodes and c corresponds to the con- 
nectance, defined as c = L/[s{s — 1)], where L is the number of 
directed links). 

Furthermore, AUesina et al.^"' generalized the semicircle law under 
the condition that A is sparse (i.e., as s ^ sc, sc k, where A: is a 
constant), and they proposed the semi-superellipse distribution: 



Fi{A = x)=P{x) = 



4(2r)Y(l-hl/n)'r(l + 2/n)-'' 



(1) 



where the parameters r and n need to be estimated using different 
equations (see^^ for details). 

Wigner's semicircle law can be considered as a particular case of 
the semi-superellipse law when n = 2. In this case in particular, the 
semi-superellipse distribution is equivalent to the semicircle distri- 
bution: 



Pr(A = x)=P(x) 



2\/(2r) 



(2r)'7t 



However, when £(.4;,) ^0, we cannot directly consider the semi- 
superellipse distribution (i.e.. Equation (1)) for estimating the real 
eigenvalue density of A for bipartite networks. In particular, a pre- 
vious study^'' emphasized the importance of a modification of the 
semi-superellipse law in order to consider non-negative matrices 
(i.e., those with entries greater than or equal to 0). 

For example, the semi-superellipse law cannot explain the tend- 
ency in non-negative matrices for the dominant eigenvalue Xx to 
increase faster than the second eigenvalue A2 in Erdos-Renyi (ER), 
scale-free, clustering random graphs''^^"'''. In particular, the spectral 
gap, defined as A1//I2) grows with increasing network complexity (i.e., 
increasing c), indicating the detachment of Ai from the rest of the 
eigenvalues. In this context, a previous study"'^ considered ~Ai sepa- 
rately from a semi-supereUipse distribution (see Equation (S23) in 



Ref 22) when estimating the dominant eigenvalue of a unipartite 
network (i.e., non-negative matrix). 

We here demonstrate that the detachment of ~Ai is also observed in 
bipartite networks. A bipartite network or graph contains s nodes 
(species) that are classified into two disjoint sets, A (animals in pol- 
lination or seed dispersal networks) and P (plants), and has L/2 
undirected links that are drawn between nodes in set A and nodes 
in set P only (i.e., there are no links between nodes belonging to the 
same set)'. Binary networks can be represented as an adjacency 
matrix A-, in which A\j = 1 if / and j are connected and 0 otherwise. 
For quantitative or weighted networks, Aij has a positive value e (0, 
°c) if species i has a mutualistic interaction with species j. The largest 
(real) eigenvalue of A (i.e., spectral radius) is also known as the 
dominant eigenvalue A]. As an example, we consider a bipartite ER 
random network that consists of P = A = 5/2 nodes and L/2 edges, in 
which L/2 edges are drawn between randomly selected plant and 
animal nodes, avoiding duplicate edges. As shown in Figure 1, the 
spectral gap Xil~A2 increases when the average species degree (fc) = LI 
(2s) increases. 

In contrast to the case of unipartite networks, bipartite networks 
are well-known to have symmetrical eigenvalue distributions. In 
particular, bipartite networks have no triangles (i.e., 3-node cliques) 
because of its definition. Thus, the traces of odd powers of A equal 
zero: tr(^^) = 0, where z= 3,5,7, .... Further, the odd raw moment 
(.1^ of the eigenvalue distribution of A equals zero because 
f.i^ = tT{A^). Therefore, the eigenvalues obtained from both binary 
and quantitative bipartite networks follow a symmetrical distri- 
bution with zero mean. 

We highlight two features of eigenvalue distributions obtained 
from bipartite networks: 1) the detachment of from the bulk of 
the eigenvalues and 2) its symmetrical distribution with zero mean. 
These properties suggest that the smallest eigenvalue — /i is also 
detached from the bulk of the eigenvalues. 

Thus, we assume that the eigenvalue distribution of A are 
described by a semi-supereUipse plus two detached eigenvalues: the 
largest eigenvalue ?.i and the smallest eigenvalue A numerical 
simulation using bipartite ER random graphs has demonstrated that 
these assumptions (i.e., features 1) and 2) explained above) in eigen- 
value distributions are suitable. As shown in Figure 2, the density of 
eigenvalues follows a semi-superellipse plus twins distribution when 
networks are sparse. When the connectance increases, the density of 
eigenvalues converges to a semicircle plus twins distribution, a par- 
ticular case of the semi-superellipse plus twins law when n = 2. 

We confirm that these assumptions are also satisfied in empirical 
mutualistic networks: all eigenvalues fall in a symmetric distribution, 
while the largest and smallest ones are found far from the mean 
eigenvalue. (Figure 3; see also Supplementary Figs. SI and S2). 

With consideration for the two features of bipartite networks, 
evaluated from numerical simulations, the previous methods^^ can 
be easily extended to bipartite networks. In particular, we propose 
two methods for estimating the dominant eigenvalue of bipartite 
networks: the Semi-superellipse Plus Twins and Semicircle Plus 
Twins methods, respectively. 

We here only briefly introduce our estimation methods. Detailed 
derivations of the following equations and a description of the sym- 
bols are given in the Supplementary Information. 

As do previous methods^'^, the semi-superellipse plus twins 
method also assumes that the bulk of the eigenvalues follow a 
semi-superellipse distribution. Since bipartite networks have a sym- 
metrical relationship between the largest and smallest eigenvalues 
(i.e., ).i and — /li), the even z-th trace is described as 



tr(^') = (5-2) x/(^ + 2a^, 



(3) 



where ^i^ denotes the z-th raw moment of a semi-superellipse distri- 
bution (Equation (1)). This equation is an extension of the previous 
method (i.e.. Equation (S23) in^"") considering the existence of twin 
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eigenvalues (i.e., li and — li). In addition, the odd z-th traces equals 
zero because of the nature of bipartite networks. 

According to the previous method"^^, we thus need to solve the 
following simultaneous equation to estimate the dominant eigen- 
value Ai. 



where 



tr[A^) = {s-2) X fi2{r,n)+2Al 
tr{A*) = {s-2) X ^^{r,n) +2a*, 
tr{A'') = {s-2)xfi^{r,n)+2Al 



(4) 



3r(i)r(!) 



(5) 



^,^ir,n)= f_:^^x^P{x)dx='-^^r^ 



In this study, we estimate Xi in addition to r and n using the Nelder- 
Mead method because Equation (4) may not be explicitly solvable. 

When n = 2, the semi-superellipse plus twins method is equivalent 
to the semicircle plus twins method. In this method, we thus assume 
that the bulk of the eigenvalues fall in a semicircle distribution. In 
particular, we need to solve the following simultaneous equation to 
estimate the dominant eigenvalue li. 

(s-2)^14-|-2Aj, 



tr(^^) 
tr(^*) 



(s-2)/i2 + 2Ai 



where jij and /I4 denote the 2nd and 4th raw moments of a semicircle 
distribution (Equation 2), respectively. 

In semicircle distributions, (12 = 1^ and 1.I4 = 2r*; thus, the solution 
of Equation (6) is as follows: 



triA") tT{Ay 2tr{A' 



(when 5>0). 



(7) 



Using the this equation, we can discuss the relationship between the 
network parameters and dominant eigenvalue because the traces of a 
matrix correspond to the network parameters as follows. 

In the binary case, tr(^A^) =L, and tr(yl'') =2H(;— L-I-C4, where 
L— _j ki. Here, /c, corresponds to the number of neighbors of 
node ! (i.e., node degree), if<: = fc? denotes the sum of squared 

node degrees, reflecting the heterogeneity of the node degrees, and C4 
is the number of four-link cycles (i.e., squares) and indicates the 
extent of clustering in a bipartite network. 

Taken together. Equation (7) can be rewritten as follows: 



C4 21 



(8) 



On the other hand, in the quantitative case, tr(^^) = W2 and 
tr(^'') =2H,- W4-|-M'4C4, where H^ = q]- Here, q, = 

corresponds to the strength of /, where n{i) is the set of 
neighbors of node Furthermore, ^2 = j ^ m'^j and 
W4= _j -1 ^'^^''^ ^"j denotes link weight of between 
nodes i and j. The average of the product of link weights in four-link 
cycles is indicated by 1^4. 

In this case. Equation (7) can be rewritten as follows: 



Ha 



Wi W4C4 

T'^ 2 



{W2) 2W2 

s 



(9) 



Comparison with empirical mutualistic networks. To evaluate the 
validity of our estimation methods, we compared the predicted and 
observed dominant eigenvalues for 40 mutualistic ecological 
networks (see Methods). We considered both binary and 
quantitative networks. In this study, quantitative networks are 
defined as preference matrices generated using the method of (see 
Methods). 
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Figure 2 | Eigenvalue distributions of bipartite random graphs. The cases 
that the average degree (k) = 1, 5, and 10 are shown. 

Similarly to a previous study''^, we also considered a simple estima- 
tion method by Chung et al.''^, in which the dominant eigenvalue in a 
binary network is described as Ai = {k^)/{k) (i.e., the average of the 
squared degrees divided by the average degree). Chung's approxi- 
mation is satisfied when the minimum degree is not too small com- 



pared to the mean degree; in addition, this approximation method 
has been independently obtained by Nadakuditi and Newman^'* 
using free probability theory. Chung's approximation method can 
be extended to quantitative cases: Ai = {cf^liq) (i.e., the average of the 
squared strengths divided by the average strength). 

Figures 3 and 4 (see also Supplementary Figs. S 1 and S2) show that 
our methods are better than Chung's method in terms of the relative 
error between estimated and observed values. In addition, the pre- 
diction accuracy of the semi-supereUipse plus twins method is higher 
than that of the semicircle plus twins method. 

The semicircle plus twins and Chung et al. methods are simpler 
than the semi-superellipse plus twins method. In particular, the 
semi-superellipse plus twins method requires a numerical method 
in order to obtain the dominant eigenvalue because Equation (4) may 
not be explicitly solvable, although the other methods can explicitly 
estimate the dominant eigenvalue. Because of the simplicity and high 
validity of the semicircle plus twins method, we use this method 
when discussing the relationship between structural properties such 
as Hq and the dominant eigenvalue in the following sections. 

Theoretical predictions for local stability of mutualistic networks. 

We here discuss the local stability of mutualistic ecosystems using 
our estimation method. 

As briefly explained in the introduction, local stability analysis is 
often useful for evaluating the resilience of a dynamical system at 
equilibrium to perturbations. According to a previous study*, in 
particular, the local stability can be interpreted as follows. When 
an equilibrium point is stable, then the system returns to that point 
despite small perturbations. For unstable equilibrium points, small 
perturbations will move the system away from the original fixed 
point. The stability at an equilibrium point is completely determined 
using the real parts of the community matrix eigenvalues. In particu- 
lar, the equilibrium point is stable when all eigenvalues are negative. 

Equation (8) shows that the dominant eigenvalue of binary bipart- 
ite networks strongly depends on the heterogeneity of node degrees 
(i.e., Hi^. Equation (9) demonstrates that the dominant eigenvalue of 
quantitative bipartite networks is strongly affected by the heteroge- 
neities of node strengths and link weights (i.e., JJ^ and Wj, respect- 
ively). In addition to this, the number of four-link cycles C4 is also a 
determination factor. 

We investigate the effect of these structural features on local 
stability (i.e., the dominant eigenvalue). In particular, null model 
analyses are useful for evaluating the effect of structural features on 
the local stability (see Methods for details). We consider four null 
models (see Methods and Supplementary source code for details), all 
randomized networks generated from an empirical quantitative net- 
work (i.e., M). In particular. Null model 1 (NMl) is used to invest- 
igate the impact of the allocation to link weights on the networks. 
Null model 2 (NM2) and Null model 4 (NM4) evaluate the effect of 
the heterogeneity of link weights W2 and the impact of the hetero- 
geneity of node degrees H^-, respectively. Null model 3 (NM3) is 
useful for discussing the influence of structural properties, excluding 
W2 and H^. 

We can evaluate the statistical significance of structural features 
using the null models. In particular, the probability (i.e., p-value) of 
observing a network feature from a null model ensemble that is 
greater in value than that of the empirical network property is useful 
for characterizing the statistical significance of network structural 
properties'". A p-value of <0.05 indicates the network measure of 
an empirical network is larger than that of null models, suggesting 
that structural properties, collapsed because of randomization, con- 
tribute to increase that measure. On the other hand, a p-value of 
>0.95 implies that such properties contribute to decrease that mea- 
sure. All p-values of network measures are calculated from 1000 null 
model networks. 
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Heterogeneity of link weights and node degrees determines local 
stability. The semicircle plus twin method (i.e.. Equation (9)) 
predicts that Hq, W2, and C4 determine the local stability when s 
and L are constant. In particular, a linear relationship between 
■\/Hci and ?Li is expected if W2 and C4 are constant. 

This hypothesis can be tested using NMl, in which network topo- 
logy and link weight distribution are preserved. As expected, we found 
approximately linear relationships between {/H^ and li in the null 
model networks (Figure 5 and Supplementary Fig. S3; the average 
coefficient of determination of 40 empirical networks is 0.84). 

However, the randomization of link weights hardly affects the 
local stability. In 36 (90%) samples, a statistical difference between 
the dominant eigenvalues of the empirical networks and NMl could 
not be concluded (i.e., 0.05 < p < 0.95). This may be because this 
randomization hardly influences Hq. In fact, we found no statistical 
difference for if^ between empirical and nuU model networks in 36 
(90%) samples (i.e., 0.05 <p< 0.95). (Supplementary Table S2). 

We expect that the limited effect of randomized link weights is 
because the link weight distributions of the empirical networks are 
identical to the null model networks. Thus, under the condition that 
network topology is preserved, we evaluated the effect of link weight 
distributions on local stability using NM2, in which heterogeneity of 
weight distributions is controlled by the parameter a of log-normal 
weight distributions (see Methods for details). In particular, it is 
expected that Ai of the empirical networks will be significantly larger 
than that of the null model networks (i.e.,p < 0.05) as a result of the 



decrease of and W2 caused by the decrease of a. As expected, a 
lower p-value for the empirical Ai in addition to and W2 was 
observed when a is smaller (Figure 6 and Supplementary Fig. S4). 

Numerical simulations using NM3 and NM4 suggest the limited 
effect of C4 on }.i in the empirical networks. In particular, linear 
relationships between \/~tiq and li were observed (see 
Supplementary Figs. S5 and S6, respectively) although C4 was not 
conserved because of the randomization of network topology. 

Furthermore, a comparison between the dominant eigenvalues of 
empirical networks and NM3 (see Supplementary Fig. S7) indicates 
the limited effect of other topological features, excluding degree dis- 
tributions. We found no statistical difference of the dominant eigen- 
value in 36 (90%) empirical networks (i.e., 0.05 < p < 0.95). The 
limited effect of the randomization of network topology (under the 
condition that degree distributions are preserved) on the local 
stability is also because this randomization hardly influences H^. 
No statistical difference of between empirical and null model 
networks was observed in 37 (92%) samples (i.e., 0.05 < p < 0.95) 
(see Supplementary Table S3). 

On the other hand, an evaluation using NM4 demonstrated a 
statistical significance of degree distributions on local stability. In 
26 (65%) samples, we found the Xi of the empirical networks was 
significantly larger than that of the nuU model networks (i.e., p < 
0.05). This is because of nuU model networks decreases (p < 0.05) 
because there is no conservation of degree distributions in empirical 
networks (see Supplementary Table S4). 
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(right) are shown. 



Taken together, these results imply that the heterogeneity of link 
weights and node degrees (i.e., link weight and degree distributions) 
mainly determine the local stability of mutualistic ecosystems. 

Local stability is determined by network heterogeneity rather than 
nestedness. Nestedness (in terms of a topological or structural fea- 
ture) is explained in binary and quantitative networks as follows. In a 
nested binary network, interactions are organized such that 
specialists (for example, pollinators that visit a few (specific) 
plants) interact with subsets of the species with whom generalists 
(for example, pollinators that visit a number of plants) interact. In a 
nested quantitative network, interactions are organized such that 
specialists interact with species more weakly than generalists 
interact with that species. 

Although the importance of nestedness in ecosystem stability has 
been reported' our methods and the above numerical simulations 
suggest a limited effect of topological nestedness on the local stability. 
Rather, the local stability is determined by Hj- in the binary case and 
ifq in the quantitative case. We test this theoretical prediction here. 

The effect of nestedness on the local stability can be evaluated 
using NM3 and NM4. To measure topological nestedness, we used 
NODF"' and CMNB' for binary networks, and WNODF^" and 
WINE'' for quantitative networks (see Methods for details). 

Although Staniczenko et al.' has demonstrated this using a differ- 
ent approach, we also found that nestedness is not a significant 
structural pattern in quantitative networks (see Supplementary 
Table S3). In particular, the quantitative nestedness of 35 (88%) 
empirical networks is equivalent to that of the null model networks 
(i.e., 0.05 < p < 0.95). 

To investigate the relationship between the dominant eigenvalue 
and network measures, we performed a correlation analysis. AU of the 



nestedness measures showed a weak correlation with the dominant 
eigenvalues obtained from null model networks, generated using 
NM3 and NM4, while the strength heterogeneity showed a strong 
correlation (Figure 7; see also Supplementary Tables S3 and S4). 

However, this result does not indicate that there is no correlation 
between nestedness and local stability (i.e., Ai). In particular, Jonhson 
et al.'^ showed that the degree distribution (i.e., ifj-) mainly deter- 
mines nestedness (see also Supplementary Table S5). In addition to 
this, a perfectly nested graph shows the largest eigenvalue when 
considering a network consisting of a given numbers of nodes and 
links'. These results suggest a positive correlation between nested- 
ness and Xi. We evaluated this prediction using NM4 in binary case. 
As expected, statistically significant positive correlations between 
nestedness and ).x were observed (see also Supplementary Table 
S5). However, it is expected that is determined by Hj- rather than 
nestedness because degree-degree correlation (assortativity) in addi- 
tion to also affects nestedness'^ In particular, correlations of X\ 
with nestedness are expected to be weaker than those with H^- because 
of variability in assortativity. In fact, such a tendency was confirmed 
in binary networks (Supplementary Table S5). Note that we only 
considered NM4 in this case because we wanted to evaluate which 
factors or binary nestedness) mainly determine X^. In this con- 
text, NM3 is not useful because is trivially constant because of the 
preservation of degree distributions. 

This theoretical framework mainly explains the case of binary 
networks; however, it may also be applicable to weighted networks. 
Using NM4, similarly, in 28 (70%) empirical networks, statistically 
significant positive correlations between X^ and nestedness 
(WNODF) were also observed while the correlation coefficients 
between X^ and WNODF were smaller than those between X^ and 
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Figure 5 | Linear relationship between ^/H^ and Ai in NMl. The red symbols indicate empirical networks. The dashed blues line corresponds linear 
regressions between \/H<^ and /li. R' denotes the coefficient of determination for a linear regression model. 



(Supplementary Table S6). A similar tendency was also observed 
using NM3 (i.e., when is changed while the degree distribution is 
preserved); in particular, the positive correlations were concluded in 
35 (88%) networks (Supplementary Table S7). In a few empirical 
networks, any correlation between li and WNODF was not con- 
cluded. This may be because of link weights. When directly applying 
the theoretical frameworks'" '^ to weighted networks, several assump- 
tions are required (e.g., link weights are drawn from a homogenous 
distribution such as a normal distribution). However, in empirical 
networks, such assumptions may be not satisfied. In this case, a 
weaker correlation between li and nestedness (WNODF) is expected 
to be observed when heterogeneity in a link weight distribution is 
higher. In fact, the link weight heterogeneity (i.e., standard deviation 
of Wij) shows a negative correlation with the correlation coefficient 
between and WNODF in both NM4 (Spearman's rank correlation 
coefficient p = —0.67 andp = 3.7 X 10"''; see also Supplementary 
Table S6) and NM3 (p = -0.71 and p = XIQ-''; see also 
Supplementary Table S7). 

These results imply that the heterogeneity is the primary factor for 
predicting the local stability of the empirical networks and that the 
topological nestedness is a secondary factor. This result may be 
related to the conclusion obtained from a previous study" that 
showed that the number of mutualistic partners of a species is a much 
better predictor of individual species survival and nestedness is a 
secondary covariate. 

Discussion 

We proposed theoretical methods for estimating the dominant 
eigenvalue of both binary and quantitative bipartite networks. 



Although the methods are obtained as natural extensions of the 
previous methods^^, we obtained an interesting result. In particular, 
we revealed a mathematical architecture in the relationship between 
ecological mutualistic networks and local stability. Our methods and 
numerical simulations clearly showed the local stability is deter- 
mined by the heterogeneities of node degrees and link weights rather 
than topological nestedness such as NODF and WNODF. This study 
is consistent with the conclusion obtained from the previous stud- 
jgg6,i4,i5. particular, it provides more conceptual (theoretical) evid- 
ence for the limited impact of nestedness in mutualistic ecosystems. 

However, this conclusion is limited to the context of topological 
nestedness such as NODF and WNODF. A previous study* has cast 
doubt on the importance of NODF and WNODF for measuring 
nestedness; in particular, it has argued that nestedness is strongly 
related to the dominant eigenvalue (i.e., local stability) according to 
the mathematical fact that the perfectly nested graph shows the 
dominant eigenvalue. That is, it remains possible that topological 
nestedness such as NODF and WNODF cannot well capture this 
mathematical feature. These facts highlight the need for a more 
suitable definition of nestedness although Staniczenko et al.' have 
proposed to directly use the dominant eigenvalue when evaluating 
nestedness. 

In addition, more careful examinations may be required because a 
number of factors influence ecosystem stability. One remarkable 
example of this is modularity. Modularity is observed in mutualistic 
networks'*^, and it is believed to decrease their persistence". 
However, we believe there is a limited effect of modularity on local 
stability in empirical networks because the previous study" reported 
a weak significance of modularity in mutualistic networks, and our 
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Figure 8 | Schematic diagram of the null models. Null model networks, generated using Null Models 1-4, from an empirical ecological network are 
shown. 



numerical simulations using null models demonstrated that the local 
stability does not change, even if the modularity- related parameter 
(i.e., C4) is not preserved. In this manner, we carefully evaluated the 
effect of structural properties and dominant eigenvalue as much as 
possible, using several types of null models; however, it remains 
possible that other hidden structural properties primarily affect the 
dominant eigenvalue. Our analysis has such limitations, as do many 
other works on network analyses. 

The definition of ecosystem stability is controvertible. Our finding 
seems to be contradict a conclusion" that higher diversity (i.e., s) and 
connectance (i.e., L/[s(s — 1)]) promote the persistence and resili- 
ence of mutualistic communities. In particular. Equations (8) and (9) 
suggest that higher network complexity (i.e., diversity X connectance 
= L/s) leads to lower local stability, similarly to May's stability cri- 
terion'. However, this inconsistency between our conclusion and the 
previous conclusion is because the ecosystem stability definitions of 
persistence and resilience are different to the definition of local 
stability. In particular, the persistence indicates the proportion of 
persisting species once equilibrium has been reached, and the resi- 
lience represents the speed at which the community returns to the 
equilibrium after a perturbation". For example, locally stable eco- 
systems can be persistent; however, locally unstable ecosystems may 
still display persistence because of the existence of alternative stable 
states". In this case, trivial local stability may be observed because the 
intra-species interactions of the community matrix are large enough 
to compensate for the potential destabilizing effect of heterogeneity. 
A future challenge is to combine local stability with persistence and 
resilience in order to analyze the overall robustness of mutualistic 
ecological communities. 

Equations (8) and (9) also imply that the heterogeneity of network 
structure and interaction strength decreases the local stability. This 



finding may answer why a weak significance of heterogeneous degree 
distributions is observed in empirical mutualistic networks'^. That is, 
mutualistic ecosystems may avoid such a heterogeneous community 
structure in order to maintain or increase local stability. In addition 
to this, the avoidance of interaction strength heterogeneity may be 
related to a biologically feasible assumption that interaction strength 
decreases with increasing number of interacting species (i.e. node 
degree)^. In such a case, interaction strength homogeneity may 
remain despite the increase of interspecific interactions. This may 
be also a strategy for increasing the local stability of ecosystems. 

These findings emphasize the importance of heterogeneity of 
mutualistic networks in ecosystem stability, and they enhance our 
understanding of structure-stability relationships. 

The spectral radius is linked to local stability and other dynamical 
functions in wide-ranging networked systems'^. The proposed frame- 
work may be useful for the theoretical analysis of a wide variety of 
systems. 

Methods 

Empirical mutualistic networks. We collected 40 empirical mutualistic networks 
from the literature and databases (see Supplementary Table SI). In particular, 22 
networks were collected from the Interaction Web DataBase (IWDB) (www.nceas. 
ucsb.edu/interactionweb), 16 networks and two networks were obtained from Refs. 
36 and 37, respectively. 

Quantitative networks. Several types of community matrix A4 (i.e., quantitative 
networks) have been proposed. For example, there is a definition in which link 
weights (i.e., number of visits) of the mutualistic network act on behalf of the value of 
the community matrix elements, based on the arguments of the strong positive 
relationship between the interaction frequency and strength"*^"^'. However, this 
definition does not consider species abundances. A weak per capita interaction 
strength of species j upon species i can stiU result in a sizable interaction strength at 
the population level if species i is abundant. Thus, in this study, we considered another 
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definition: preference matrix A4, attained from the original mutualistic visiting 
network after adjusting for uneven species abundances under the mass action 
hypothesis (see^ for details). 

Null models. Randomized graph ensembles, in which a part of structural properties 
in an empirical network are preserved, are often used as nuU models. 

Mutualistic networks are strongly influenced by species abundance, trait matching, 
spatio-temporal variation of individuals and species, phylogenetic relatedness, and 
sampling efforts'*^ '*^. Such effects determine not only the existence of interactions but 
also the frequency of these interactions among species. Thus, we proposed four null 
models by considering two aspects in randomization: (1) randomization of interac- 
tions among species (i.e., switching or rewiring links) and (2) randomization of 
interaction strengths (i.e., reshuffling or randomly reassigning link weights). Figure 8 
illustrates a schematic diagram of the null models. 

In Null model 1 (NMl), the binary structure (i.e., topology) of the empirical net- 
works is preserved; however, the link weights are switched between two randomly 
selected links until the switching of all link weights is completed. 

In NuU model 2 (NM2), as in NMl, the topology of the empirical networks is 
preserved; however, the link weights are drawn from a probability distribution under 
the constraint that the total of link weights is equivalent between the empirical and 
null model networks. In this study, inspired by previous studies*^^-^^, we assume that 
the link weights in a community matrix are drawn from a log-normal distribution 
(i.e.. In W-^A/^f/^o"^)). To consider the above constraint condition, we used an 

algorithm proposed in Ref. 49. Since E[W] - e^'^-^"' and Var[W] - (^e^' - 1 j e^^'^"' - 

(e'^^ — l)(E[W])^, <7 indicates the variance of log-normal distributions if E[TV] is 
constant. Thus, we can control the variance of log-normal distribution Var[VV^], 
which is related to the heterogeneity of link weights W2, using o only. 

In NuU model 3 (NM3), links among species are rewired randomly while con- 
serving the degree distribution of empirical networks. In particular, we used an edge 
switching algorithm"*" "*^ that generates a random network by switching two randomly 
selected links until the switching of aU links is completed. Note that the link weights 
are also implicitly switched according to its link switching. 

SimUarly to NM3, NuU model 4 (NM4) also considers the randomization of net- 
work topology; however, the degree distributions of empirical networks are not 
conserved. This nuU model corresponds to the bipartite ER random graph model. In 
particular, the nuU model networks were generated to maintain the numbers of 
plants, animals, and the weighted links in empirical networks. 

Nestedness. We computed four nestedness measures for binary networks and 
quantitative networks: 1) Nestedness based on Overlap and Decreasing FiU 
(NODF)^^, a popular nestedness measure for binary mutualistic networks; 2) 
Nestedness based on Common Neighbors (CMNB), defmed by BastoUa et al.^; 3) 
Weighted NODE (WNODF)-'°, a straightforward extension of NODE to quantitative 
networks; and 4) Weighted-Interaction Nestedness Estimator (WINE)-^' that uses the 
weighted Manhattan distance in order to calculate quantitative nestedness. 

AU these nestedness measures except CMNB were calculated in R version 3.0.2 (R- 
project.org) using the nested function in package bipartite version 2.02. 
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